Ground-state properties of the Rokhsar— Kivelson dimer model 

on the triangular lattice 
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We explicitly show that the Rokhsar-Kivelson dimer model on the triangular lattice is a liquid 
with topological order. Using the Pfaffian technique, we prove that the difference in local properties 
between the two topologically degenerate ground states on the cylinders and on the tori decreases 
exponentially with the system size. We compute the relevant correlation length and show that it 
equals the correlation length of the vison operator. 
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I. RVB LIQUID AND THE 
ROKHSAR-KIVELSON DIMER MODEL 

Resonating valence bond (RVB) spin liquid in two di- 
mensions is a remarkable theoretical concept which pre- 
dicts very unusual low-temperature properties of spin-1/2 
systems jlj. Unlike conventional magnetically-ordered 
ground states with spin-1 excitations, the RVB ground 
state has no long-range order of any local order pa- 
rameter and possesses elementary excitations with spin 
1/2|Q. If such a system is doped with mobile holes, 
the fractionalization of spin excitations translates into 
the effect of spin-charge separation: this scenario has 
been widely explored in the context of high-temperature 
superconductivity)!, ||, (§• While a rigorous verification 
of the RVB liquid phase in a realistic spin system is usu- 
ally very difficult due to the strongly-correlated nature 
of the state, many specially designed systems have con- 
firmed RVB liquid properties § |, g § g, 0. 

In understanding generic properties of the RVB liquid 
state, one may benefit from studying dimer models which 
are closely related to the RVB spin liquids pd|, p"2[ . In the 
RVB construction, the wave function is represented as a 
sum over singlet configurations, while in the dimer mod- 
els the singlets are replaced by dimers. The difference be- 
tween the spin and the dimer systems is that dimer con- 
figurations are mutually orthogonal by definition, while 
different singlet configurations have a finite overlap |]l3|. 
The RVB spin liquid is known to have two types of ele- 
mentary excitations: spinons (spin-1/2 excitations) and 
visons (Z12 vortices) 14]]. In a dimer liquid, the spinons 
are prohibited by the dimer constraint (or, equivalently, 
are pushed infinitely high in energy) , while the visons are 
expected to be indeed the lowest excitations above the 
ground state. The RVB spin liquid must have topologi- 
cal degeneracy on domains of nontrivial topology. Such 
a decoupling into topological sectors is straightforward 
in dimer liquids where it is defined in purely geometric 
terms 0, [l^, Thus dimer models provide a conve- 
nient test ground for studying properties of RVB liquids 
related to vison excitations. 

The advantage of studying dimer models is that they 
are simpler than spin models, and are hence better un- 
derstood and more accessible to analytical methods. The 



simplest dimer model contains the pair-wise hopping 
term and the pair-wise potential term 
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where the sum is taken over the four-vertex plaquets of 
the lattice |lf|. On the square lattice such a model has 
a crystal ground state (with the crystal of dimers break- 
ing the translational symmetry of the lattice) for any ra- 
tio v/t, except for the phase-transition point (s) between 
different crystal structures [l5|, [l7j. The situation is 
different on the triangular lattice (with the sum in (|l|) 
taken over all rhombi consisting of two neighboring tri- 
angles) : there the dimers are believed to form a liquid in 
a finite range of the parameter v/^H ^f)- Out of all pos- 
sible values of v/t, one is special: v/t — 1. At this ratio 
of parameters — further called Rokhsar-Kivelson (RK) 
point — the ground state of the Hamiltonian ([!]) is known 
exactly: it is a superposition of all possible dimer configu- 
rations with equal amplitudes [unfortunately, the excited 
states are not known exactly even at the RK point]. At 
v/t > 1, this state immediately yields to the staggered 
crystal state via a first-order phase transition. However 
at v/t < 1, from the available numerical evidence it fol- 
lows that the RK state on the triangular lattice smoothly 
evolves without crystallization in a finite range of v/t (ap- 
proximately until v/t sa 0.6 . . . 0.8) |^, It is this region 
of the parameter space which contains the dimer liquid. 

The RK point is one of the representatives of the liq- 
uid phase and is most accessible for analytic treatment, 
because averaging over the ground state is equivalent to 
the statistical averaging over all possible dimer configu- 
rations. Such an averaging may be performed with the 
usual Pfaffian technique which allows to establish the 
exponential decay of dimer correlations^]. In this paper, 
we employ the Pfaffian technique to explicitly compute 
the correlation length involved in the ground-state cor- 
relation functions and in the splitting between the topo- 
logical sectors on the cylinder and on the torus. We fur- 
ther demonstrate that the same correlation length also 
appears in the vison correlation function. Thus our cal- 
culations prove that the RK ground state on the triangu- 
lar lattice possesses topological order: on a topologically 
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nontrivial domain, the ground states in different topolog- 
ical sectors are degenerate and locally indistinguishable 
(in the limit of infinite system size). Note that the ex- 
ponential decay of ground-state correlation functions is 
consistent with the numerical finding of the gap in the 
excitation spectrum at the RK point (the magnitude of 
the gap is numerically estimated as 0.1t)|18, 20|. 



II. 



PFAFFIAN TECHNIQUE AND 
CORRELATION LENGTHS 



B 




In the Pfaffian technique, the statistical averaging over 
the dimer configurations is performed via introducing an 
auxiliary real fermionic variable (a Majorana fermion) on 
each lattice site[[l9|. The total number of dimer configu- 
rations may then be written as the partition function of 
these fermions 



Z = 



jPJ da% cxp ^""^ a,jA 



= Piaff(Ay), (2) 



where i and j label the lattice sites, and the fermionic 
variables ai obey the conventional rules: a^a, = — Oj-ai, 
/ dcii = 0, J CLidai = 1. The "hopping amplitudes" Ay- 
take values ±1 on nearest-neighbor sites and otherwise, 



and form an antisymmetric matrix: Aj 



The 



signs of Aij must be adjusted so that all terms in the ex- 
pansion of the exponent in (^) give positive weight. Such 
terms in the expansion of the exponent are in one-to- 
one correspondence to the dimer coverings of the lattice. 
Each of them has magnitude one, and should they all 
have equal signs, the partition function (j^) simply counts 
the total number of dimer coverings. 

As shown in Ref. |l^, a necessary requirement for the 
proper relative sign of different dimer coverings is that 
the circular product Jl^k; equals (—1) around any el- 
ementary rhombus (around any elementary even-length 




FIG. 1: Topological sectors on a cylinder and on a torus. 
Dimer coverings may be classified according to the parities of 
the number of dimers intersecting the reference lines (dashed 
lines). Dimer configurations differing by a rearrangement in 
contractible domains (shaded areas) belong to the same topo- 
logical sector and contribute with the same sign to the parti- 
tion function (Q). To change the topological sector, a circular 
permutation of dimers along a topologically nontrivial con- 
tour (shown in a zig-zag line in the case of the cylinder) is 
necessary. 



FIG. 2: One possible choice of the amplitudes Aij. This 
choice of amplitudes is periodic with a unit cell containing 
two lattice sites (marked by a dashed ellipse). The arrow 
directions correspond to the signs of Ay: Aij equals 1 if the 
arrow points from i to j, and equals —1 if it points from 
j to i. Also shown are A and B directions: parallel and 
perpendicular to the lattice lines, respectively. 



cycle, in the general case of a planar lattice). Then it 
follows that Y[ = — 1 over any contractible contour of 
even length and, as a consequence, any two dimer con- 
figurations differing by a local rearrangement of dimers 
contribute to the partition function (||) with equal signs. 
Different dimer configurations which cannot be related 
by local dimer rearrangements (i.e. by rearrangements in 
a contractible domain, see Fig. Q) may contribute with 
either equal or opposite signs: this effect will be of crucial 
importance for our calculation. 

One 
in Fig. 



possible choice of the amplitudes Ay is shown 
4 Note that in order to obey the sign rule for 
A^ we have to double the unit cell of the lattice. How- 
ever, any physical quantity computed with those Aij has 
the periodicity of the lattice, as the set of amplitudes 
A.^ translated by one lattice spacing may be returned to 
the original form by an appropriate Z 2 gauge transfor- 
mation Aij I— > WiAijWj with Wi — ±1. Now the par- 
tition function (^]) as well as correlation functions may 
be conveniently computed in the Fourier components. 
The Fourier transformation of Ay is a 2x2 matrix A(k). 
From antisymmetricity of the matrix Ay it follows that 
its eigenvalues always come in pairs ±i-E(k), with com- 
plex conjugate eigenvectors corresponding to opposite 
eigenvalues. A straightforward calculation shows that, 
for the triangular lattice, the spectrum of A(k) in the 
bulk is always gapped (never crosses zero as a function 
of k) . Explicitly, it is given by 



E(k) = 2 (cos 2 ki + cos 2 k 2 + cos 2 k 3 



1/2 



(3) 



where k\, k 2 , and k% are the projections of the vector k 
onto the three lattice directions (obeying the constraint 
k\ + k 2 + &3 = 0). We have also shifted the origin of 
the Brillouin zone to bring (^) to a symmetric form. The 
gap in the spectrum of Ay is crucial for the exponen- 
tial decay of correlation functions: as we shall see below, 
the correlation length is determined by the complex wave 
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vectors k solving the equation 2£(k) = 0||. At this stage 
of calculation, we see the difference between the triangu- 
lar and square lattices: on the square lattice, the corre- 
sponding matrix Aij leads to a gapless spectrum (and, 
consequently, to the power-law decay of correlations pjj . 
Note that although the presence or absence of the gap in 
the matrix Aij (governing the decay of correlation func- 
tions in the ground state) agrees with the presence or 
absence of the gap for excitations, there is no direct re- 
lation between the two gaps. The spectrum of is not 
the energy spectrum of quantum dimers, but only the 
spectrum of auxiliary Majorana fcrmions introduced for 
calculating the classical partition function of all dimcr 
configurations. 

Correlation functions of dimers may be expressed in 
terms of Green's functions of Majorana fermions. In- 
deed, the probability of dimers occupying a given set of 
positions (connecting pairwise points 1 and 2, 3 and 4, 
. . . , 2m — 1 and 2m) may be computed by excluding 
those points from the lattice, which is equivalent to plac- 
ing Majorana fermions at those points: 



(«12 



n2m-l,2ra) — {0,1(12 ■ 



■ a 2r . 



Z 1 I \\_dai (<zia 2 • ■ ■ »2m) exp 



/ J &iAijl 



(4) 



(up to a sign). Here the first average denotes the sta- 
tistical averaging of dimer occupation numbers riij over 
all dimer configurations, while the second average is the 
quantum-mechanical correlation function in the theory of 
Majorana fermions. This correlation function may fur- 
ther be decoupled using the Wick theorem. 
The Green's functions are obtained as 



G(i,j) 



{aidj) = 



dk J 4- 1 (k)e ik(r 



i) 



(5) 



.z. 



with the integral over k defined as averaging over the 
Brillouin zone. The decay at large distances is found from 
deforming the integration domain into complex values 
until it reaches the singularities in .<4 _1 (k). The singu- 
larities are the square-root branching points, which trans- 
lates into the asymptotics G(R) oc R^ 1 / 2 exp(— R/£ ). 
The correlation length £ is direction-dependent. In the 
direction specified by a unit vector n, the correlation 
length £ n is given by £~ = — «k n n, where the complex 
vector k n is determined from the set of equations: 



£(k„) = 0; 
V£(k n ) || n. 



(6) 
(7) 



(the meaning of the second equation is that the "group 
velocity" is directed along the vector n). Equivalcntly, 
the conditions (||), ^ for k n may be expressed as 



= maxmin Rc (— ikn) 



£(k)=0' 



(8) 



where maxmin Ro denotes minimizing the real part 
(among its positive values) over the components of k 
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FIG. 3: Correlation length and its direction dependence, (a) 
Construction of the correlation length. In the (Im k x , Im k y ) 
plane, there is a region around Imk = where there are 
no solutions to the equation E(k) — 0. The thick solid line 
denotes the boundary of this region [the axis k x is chosen 
along one of the lattice directions (A direction), and the axis 
k y is perpendicular to it (B direction); if reflected in those two 
axes, the thick solid line has the hexagonal symmetry]. Given 
the direction n, the vector Imk n is chosen on this boundary 
to maximize its projection onto n [in accordance with (0), 
(0), or, equivalently, with (Q)] The vector k n determines the 
correlation length 1 = — ink n in the n direction. For a vison 
tunneling as a wave front [Eq. (p^)] , the tunneling amplitude 
has a different correlation length given by (£n) _1 = —ink* 
with the imaginary component of the vector k*, parallel to n. 
(b) The direction dependence of the real parts of the inverse 
correlation lengths and (£n) -1 - The inverse correlation 
lengths are plotted versus the angle between the direction 
n and the lattice lines. The angle corresponds to the A 
direction, and the angle n/6 corresponds to the B direction. 



parallel to n and maximizing it over the perpendicular 
components. The vectors k are confined to the complex 
surface E(k) — 0. The construction of is illustrated 
in Fig. | 

In the particular cases of directions A and B paral- 
lel and perpendicular to the lattice lines respectively (as 
shown in Fig. |^), the correlation lengths may be found 
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analytically (from solving (||) and (R)) 




0.83 ± 0.12 in, 
1 



^ 1 = 7f ln ( 2 + ^ 



0.76 



(9) 
(10) 



(in the units of the lattice constant). The imaginary 
part in £7 implies the clamped oscillating asymptotic 
behavior (with an incommensurate wave vector) of the 
Green's function and of other relevant correlations (e.g., 
vison correlation function, see Section |V|). 

For other directions, the equations (|6|) and (0) [or, 
equivalently, (||)] may be solved numerically. In Fig. [|, 
we plot the direction dependence of the inverse coherence 
length f- 1 . 

The exponential decay of the Green's functions implies 
an exponential decay of all dimer-dimer correlations. For 
example, the pair-wise dimer correlation {n^n^i) decays 
with the correlation length £/2, since it involves a prod- 
uct of two Green's functions [| 



III. TOPOLOGICAL DEGENERACY ON 
MULTIPLY-CONNECTED DOMAINS 

If one considers dimer coverings on a domain of non- 
trivial topology (e.g., torus, cylinder, domains with holes, 
etc.), such coverings may be classified into topological 
classes |l4|, pi. Dimer coverings from different topologi- 
cal classes cannot be transformed into each other by local 
rearrangements of dimers. The topological classes are de- 
fined by the parity of the number of dimers intersecting 
a given non-contractible closed contour. On the cylinder 
there are two topological classes, on the torus there are 
four of them (Fig. ||). 

In the quantum dimer model, the Hilbert spaces 
spanned by different topological classes are not mixed 
by the Hamiltonian, and each topological class possesses 
its own ground state. At the RK point, each of these 
ground states has zero energy. In the general case of a 
topological liquid, the splitting of energies between differ- 
ent topological sectors must be exponentially suppressed 
in system size[|| [llj]. We verify this property later near 
the RK point to the lowest-order perturbation theory in 
(v/t-l). _ 

At the RK point, we compare correlation functions in 
different topological sectors. We shall see that the corre- 
lators in different sectors nearly coincide, with the split- 
ting exponentially small in system size. We start our 
analysis with comparing partition functions for different 
topological sectors. The derivation looks different for the 
cases of the torus and of the cylinder. 

Consider first the topological sectors on the torus. De- 
note the numbers of dimer coverings with even/odd in- 
tersection indices at reference contours by iV ee , N eo , N oe , 



and N DO , where N Ca . e corresponds to the sector with the 
parity e x of intersection with the line parallel to y-axis, 
and the parity e y is of intersection with the line parallel 
to x-axis. In the Pfaffian technique, different topologi- 
cal sectors may be accessed by imposing either periodic 
or antiperiodic boundary conditions across the reference 
contours (in the x and y directions) JT^, [l6, 22 . The cor- 
responding partition functions Za^o (cF x ,a y = ±) are 
expressed in terms of N ea . e as 



Z++ = N ee + N eo + N oe - N QO ; 
Z+- = N ee ~ N eo + N oe + N ao ; 
Z-+ = N ee + N eo - N oe + N OD ; 
Z— = -N ee + N eo + N oe + N OD . 



(11) 



(up to the simultaneous change of signs of a x and/or a y 
in all the four expressions above: the particular choice of 
signs of a x and a y depends on the gauge choice for Aij\ 
the overall sign of Z axlJy is ignored, as usual). These ex- 
pressions are a particular case of a more general formula 
for surfaces of arbitrary genus |23| . We now establish that 
the splitting in is exponentially small in the system 
size, then the same exponential smallness will follow for 
the splitting of N est . ev . Consider, for example, the parti- 
tion functions Z ++ and Z^ Each of them is given by 

the sum over the Brillouin zone 



cxp Vln£(k), 



B.Z. 



(12) 



where the lattice of points k is determined by the bound- 
ary conditions. In the partition functions Z ++ and Z^ , 

the lattices of k are shifted by half a period in the y- 
direction. We can rewrite the relative difference A in 
those partition functions as a sum over trajectories with 
odd windings in the y-direction (taking the Fourier trans- 
form of In -E(k) and using the Poisson summation for- 
mula) : 



A = 



Z-i — L — Z4 



Z 



5>£(k)-]r'ln£(k') 



B.Z. 



B.Z. 



(13) 



R 



where f(R) is the Fourier transform of lni?(k), and the 
sum is taken over all closed trajectories with odd wind- 
ings in the y-direction (see Fig. ^a). The function /(R) 
decays exponentially with distance, with the same cor- 
relation lengths £ n as the Majorana Green's functions 
(the correlation lengths determined by (j|) and (0)). The 
physical interpretation of those exponents are the qua- 
siclassical amplitudes of the vison tunneling around the 
torus Q. The splitting ( |l3| ) is dominated by the largest of 
the exponents, i.e., by the optimal trajectory (the choice 
of the optimal trajectory depends on the aspect ratio of 
the torus, and, because of the anisotropy, may not nec- 
essarily be the geometrically shortest one). This result 



may formally be written as 



A < exp — min Rc 

' ' R 



R 

in 



(14) 



where £ n is given by (||) and (0) [or, equivalently by (j|)] 
with the vector n in the direction of R, and the minimiza- 
tion is performed over all displacements R returning to 
the original point with an odd winding number in the 
y-direction (Fig. Ha). 

Note that Eq. (|14|) neglects interference between dif- 
ferent trajectories in Fig. ^a. Such an interference may 
additionally suppress the splitting, hence we put the < 
sign in (U4). We can take the interference into account 



by performing not the full Fourier transformation in (13), 
but only the Fourier transformation in the y direction, 
leaving the k x component of the wave vector quantized. 
We denote this quantization of k x as k x S A, where A 
is the set of allowed values of k x [necessarily real!] de- 
pending on L x and on the boundary condition. Then 
repeating the above argument, we arrive at the estimate 
for the splitting: 



A < exp ^— min Im kL 



V \E(k)=0 



(15) 



where the minimum is taken among positive values only. 
Because of the constraint k x S A, the expression ( |l5|) is 
invariant with respect to L y L y ± L x , as expected [any 
displacement R in Fig. ^a with winding number one may 
be chosen as L y ]. 

We have obtained the two estimates ( |l4] ) and (|l5|) both 
derived in the limit L y 1 . One can verify that the esti- 
mate ( |l4| ) is stronger than (|l5|) when L x > L y 1, while 
( |l5| ) is stronger than (|l4| ) when L y ^> L x ~ 1 (quasi-one- 
dimcnsional limit). In the intermediate regime L y 3> 
L x 3> 1, both (|l4|) and (O) reduce to the same result 



A ~ exp ( — Rc 



min Re ( 



-ink 



(16) 



i£(k)=0 
llmL x k=0' 



Here n is the vector perpendicular to L x and L y = nL y . 
The minimization is performed over the wave vector k 
with its imaginary component parallel to n. Thus defined 
correlation length £* corresponds to the vison tunneling 
as a "wave front" (as a plane wave in the x direction) and 
is slightly different from £ n corresponding to the tunnel- 
ing of a "point" vison, see Fig. || (this difference is due 
to the lattice anisotropy). 

The splitting between Z ++ and Z^ translates into the 

exponentially small splitting between N eo and N OQ given 
by the same expressions (fj^)-(|l6|). Note that a vison 
tunneling in the y-direction relates the topological sec- 
tors with different parity of intersections with a contour 
running in the y-direction. Similar expressions apply for 
splitting in other pairs of topological sectors. 

The splitting of the two topological sectors on a cylin- 
der (or on a disc with one hole, which is topologically 
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(b) 



FIG. 4: The vison tunneling trajectories R contributing to 
the splitting between partition functions on the torus and 
on the cylinder, (a) The torus is constructed by identify- 
ing points on the plane differing by multiples of the two basis 
vectors and L y . Equivalently, it may be described by iden- 
tifying opposite sides of the unit cell (shaded parallelogram). 
The vectors R contributing to the splitting between Z++ and 
Z+- (or, equivalently, between N eo and N 00 ) [Eqs. (113) and 
(14)] belong to the lattice generated by L, x and L H and have 
odd winding numbers in the y direction, (b) The tunneling 
trajectories R determining the splitting between the topolog- 
ical sectors on the cylinder as given by Eq. (|l^) . The cylinder 
is defined by identifying points differing by a multiple of L ffl . 



equivalent to a cylinder) is derived in a way slightly dif- 
ferent from our analysis of the torus. Consider a semi- 
inhnitc cylinder (or a plane with one hole of arbitrary size 
and shape). On such a system, we consider two differ- 
ent boundary conditions for Majorana fermions: periodic 
and antiperiodic as we go around the cylinder (around 
the hole, respectively). To those boundary conditions, 
there correspond the two matrices and A- , •* . Below 
we prove that one of those matrices has a zero eigenvalue 
(with the eigenvector localized near the boundary), and 
the other does not (in the general case). 

The proof consists of four steps. First, we consider 
a semi-inhnite cylinder of a particular geometry, with 
the straight boundary parallel to one of the lattice lines 
(Fig. ||a). For such a cylinder, the spectrum can be ex- 
actly computed, and the existence of the zero mode in one 
of the two topological sectors is easily verified (for this 
particular geometry, the zero mode is strictly localized at 
the boundary row of sites) . Second, we note that the bulk 
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FIG. 5: (a) A particular geometry of the semi-infinite cylin- 
der with the straight edge in the A direction. For this cylinder, 
the zero mode of the Ay matrix (for one of the two boundary 
conditions) may be explicitly found. Arrows show identifi- 
cation of the two edges of the semi-infinite stripe, (b) The 
same cylinder as in (a), with a hole sufficiently far from the 
edge. Two different boundary conditions are possible across 
the dashed line, (c) A cylinder of arbitrary geometry with 
a hole far from the edge. Two different boundary conditions 
are possible across the dashed line. 

spectrum of Ay is gapped, and therefore there are only a 
finite number of states at small energies, and they are all 
localized near the boundaries. Since all eigenvalues of the 
antisymmetric matrix come in pairs HE, the property of 
having a zero mode depends on whether there are even 
or odd number of subgap states. This parity can not be 
changed by any local deformation of the antisymmetric 
matrix Ay and thus the zero mode (or its absence) is 
topologically stable (note that this argument resembles 
the proof of the topological stability of the zero mode 
in vortices in p-wave superconductors [p4|, |25f|). Third, 
once the theorem proven for a particular geometry of the 
cylinder, we make a hole in this cylinder sufficiently far 
from the boundary (Fig. ^b) . Then we can impose either 
periodic or antiperiodic boundary conditions across a line 
connecting the hole to the edge of the cylinder (dashed 
line in the figure). Switching the boundary conditions 
toggles on and off the zero mode at the cylinder edge. 
However, since the total parity of the subgap states (in- 
cluding both states at the cylinder edge and at the hole) 
is conserved, switching the boundary conditions across 
the dashed line also toggles the zero mode at the hole 
boundary. Thus we extend our theorem to the case of 
a hole in a plane (a sufficiently large cylinder may be 
regarded as a plane, from the point of view of states 
localized near the hole). Fourth, we repeat the previous 
argument for a cylinder of arbitrary geometry with a hole 
in it (Fig. ||c), and thus establish the presence of the zero 
mode at the cylinder edge for one of the two boundary 
conditions. This completes the proof. 

It is the zero modes which determine the splitting be- 
tween the topological sectors on a cylinder. The partition 
functions of the Majorana fermions on the cylinder with 
periodic and antiperiodic boundary conditions are 

Z+ = N e + N 0l Z_ = N e - N (17) 

(or vise versa, depending on the choice of the gauge for 
Aij), where N e and N a are the numbers of dimer cover- 
ings in the even and odd sectors. The splitting between 
N e and N a is given by A = Z-/Z+ and is determined 



by the splitting of the two zero modes (contributing to 
Z-). This splitting is in turn determined by the expo- 
nential tails of the zero modes away from the boundary: 
the exponential dependence of the splitting on the cylin- 
der length has the same correlation length as the decay 
of the zero modes. Similarly to the case of the torus, 
this exponential decay may be estimated in two ways. 
In the first approach, the zero mode is constructed as 
a linear combination of Green's functions localized near 
the boundary. Explicitly, we consider the zero mode 
on the "unfolded" cylinder (Fig. ||b) and act on it with 
the matrix Ay defined on the whole plane. The resulting 
wavefunciton contains positions and amplitudes of the 
sources for the Green's functions to reproduce the orig- 
inal zero mode 'J'o- This representation proves that the 
decay of zero modes (and hence the splitting A) obeys 
the estimate ([Til ) with the minimization performed over 
all possible vectors R connecting the opposite edges of 




FIG. 6: The vison tunneling trajectories determining the 
splitting between the topological sectors in the plane geom- 
etry (a disc with a hole), (a) In the general situation, the 
splitting is determined by Eq. ( |14| ) with the minimization over 
trajectories R connecting the hole to the external boundary. 
The splitting of the local dimer correlation functions (^) is 
determined by the two tunneling events: vison tunneling from 
the inner boundary to the location of the involved dimers 
(.Rin) and then further to the outer boundary (Rout)- The 
two possible situations are shown: the splitting at locations 
near the optimal tunneling trajectory have the same exponen- 
tial dependence as the direct tunneling from the inner to the 
outer boundary (dashed line), while at locations far from the 
optimal tunneling trajectory, the splitting is smaller (dotted 
line), (b) In the special case of polygon boundaries with their 
segments parallel to lattice directions (A directions), only tra- 
jectories connecting corners of the boundaries should be taken 
into account in (jlj) and (^o|), since vison tunneling from the 
boundaries is suppressed. 
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the cylinder (Fig. |^b) [see also the analogous discussion 
below of the zero modes in the plane geometry]. Al- 
ternatively, we may represent the zero mode as a linear 
combination of plane waves quantized in the ^-direction 
and decaying in the y direction, which leads to Eq. (15). 
The relation between those two estimates is the same as 
for the case of the torus, with the same expression Jl6| ) 
valid in the regime L x 3> L y 3> 1. 

We need to remark, however, that in the case of the 
cylinder, the formulas (|l4|)-(|l6|) do not properly take 
into account interference of tunneling trajectories start- 
ing or ending at neighboring points of the boundary. As 
a consequence, the expressions (^4|)-(|l^) give good (most 
likely exact) asymptotic estimates in the case of ragged 
boundaries, but may strongly overestimate splitting in 
the case of a regular edge. For example, at the straight 
edge shown in Fig. [|b, the zero mode is exactly localized 
at the boundary layer of sites, and therefore the splitting 
of the topological sectors on a cylinder with such an edge 
is exactly zero, in contrast with jhl|)-([l^). 

Our discussion of the cylinder may be directly ex- 
tended to the plane geometry (holes in the disc), with 
the only difference that there are no counterparts of 
Eqs. and (|l|) in this case. Consider specifically 

the simplest possible geometry: a disc with a hole in 
it (Fig. ^a). A reference line defining the two topologi- 
cal sectors is drawn to connect the two boundaries (the 
external boundary and the hole boundary). Similarly to 
the cylindrical geometry, the splitting of the partition 
functions Z e and Z Q is given by the splitting of the zero 
modes at the two boundaries. Our proof of the existence 
of zero modes may be easily extended to show that in 
the plane geometry, the hole encircling an odd number of 
sites hosts a zero mode for the periodic boundary condi- 
tions around the hole, while a hole with an even number 
of internal sites has a zero mode for the anti-periodic 
(vison-like, see the last section of the paper) boundary 
conditions. For a hole with an odd number of internal 
sites, the zero mode may be constructed as a linear com- 
bination of Green's functions with sources located near 
the boundary (the proof of this statement is similar to 
that in the cylindrical geometry: it follows from acting 
on the zero mode with the operator Ay defined on the 
whole plane including the hole interior). Therefore the 
zero mode decays with the same correlation lengths £ n 
as the Majorana Green's function. The same result may 
also be proven for a hole with an even number of internal 
sites: the proof is a straightforward generalization of the 
theorem about the correlation length of the vison opera- 
tor proven in the last section of the paper. Thus, for any 
type of hole, the splitting between Z e and Z a is given by 
Eq. (pT[), with the minimization performed over all tra- 
jectories R connecting the inner and the outer boundaries 
(Fig.ga). 

Note however that for certain special geometries - 
in particular, those involving straight edges stretching 
along the A directions (along the lattice lines, Fig. |^) - 
Eq. (14) may overestimate the splitting. An analysis of 



the specific case of plane geometries with only straight 
edges along the A directions reveals that tunneling of 
visons from the straight edges is suppressed (see also a 
similar remark about the zero modes at straight edges 
of cylinders), and that the splitting is given by Eq. (nj) 
with the tunneling trajectories R drawn between corners 
of the boundaries (Fig. pr 



IV. LOCAL DIMER CORRELATIONS AND 
GROUND-STATE ENERGY SPLITTING UNDER 
PERTURBATION 

Our discussion of the splitting in the partition func- 
tions may be extended to compute the exponentially 
small splitting in correlation functions of local operators. 
Consider first the simplest correlation function — the av- 
erage dimer density on a given link (ny) on a cylinder 
or on a disc with a hole. In the +/— sectors, this ex- 
pectation value may be found as the Majorana Green's 
function G y = (aiCLj). The partition functions and the 
Green's functions with periodic/antiperiodic boundary 
conditions will be denoted as Z + , Z-, and as Gf-, G~, 
respectively. As explained in the previous section, with 
one of those boundary conditions, the matrix Ay has zero 
modes at the two boundaries. Without loss of generality, 
we assume that the zero modes appear in the "— " sector. 
If we consider only one boundary (a semi-infinite cylin- 
der or only one hole), the zero mode at this boundary 
may be chosen to be real. We denote such zero modes 
by ^a{i) and ^>b{i) (a and b label the two boundaries 
near which the zero modes are localized) . In a finite sys- 
tem, these two wave functions are hybridized and split. 
The two conjugate wave functions = ^f a + i^b and 
^2 = — i^b have small non-zero eigenvalues HEq 
[the value of Eq is exponentially small in system size and 
obeys the estimates (jl4j)-([L6|)]. The difference between 
the expectation values of n y in the odd and even sectors 



Z+ - Z- 



2Z+Z. 



r^2 _ g2~(^ij ^*ij> 



2Z_ 



Z+ 



(G+- - Gj •) 



(18) 



(in the last equality we have used Z_ <C Z+). 

Further, the Green's functions may be rewritten as 
sums over eigenvectors Gy = J^n ^n{i)^n{j)/iE n . The 
main contribution to Gf, —Gj- comes from the zero modes 
and ^2 in G~. The exponentially small prefactor 
is canceled by the exponentially small denomina- 
tor Eq in G~- , and to the exponential precision we find 

{n l3 )o ~ (nij)e ~ *a(0*6(j) - *bW*a(j) (19) 

Graphically, this result may be depicted as the sum of 
the two diagrams in Fig. ^a. 

This diagrammatic approach may be extended to the 
even-odd splitting in higher-order correlation functions. 
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FIG. 7: Diagrams involved in splitting between correla- 
tion functions in different topological sectors, (a) The two 
diagrams contributing to the splitting of the average dimer 
density (riij) between the odd and even sectors on the cylin- 
der. The two edges of the cylinder are labeled a and b. The 
dashed lines represent the zero-mode wave functions ^ a and 
$(,. (b) One of the diagrams contributing to the splitting 
of the four-point correlation function {riijUki) on the cylin- 
der. The dashed lines represent the zero-mode wave functions 
ty a and ^i,; the solid line denotes the bulk Green's function 
Gjk- (c) One of the diagrams contributing to the splitting 
of the four-point correlation function (riijUki) on the torus. 
The thick solid line denotes the bulk Green's function Gu. 
The dashed line denotes the Green's function with a winding 
around the torus Gj,fc+R, where the winding vector R is one 
of those shown in Fig. 



The diagrams are drawn by pairing all Majorana fcrmions 
in the correlation function (Q) according to the Wick rule, 
except for one pair which is connected to the boundaries. 
To "internal" couplings, we associate the bulk Green's 
functions, while for the "external" lines, we associate the 
zero-mode wave functions as in (|w|) (an example of such 
a diagram is shown in Fig. 0b). 

Furthermore, the same type of diagrammatic rules may 
be derived for the case of torus, with the only difference 
that the the "external line" must in this case be closed 
and associated with the Green's function Gy+R,, where 
the vector R winds around the torus with the appropriate 
winding number (as in Fig. ^a). An example of such a 
diagram is shown in Fig.^Jc. 

In the case of cylindrical and plane geometry, the 
smallness of the correlation-function splitting follows 
from the exponential behavior of the boundary zero 
modes. Since the zero modes rapidly decay away from 
the boundaries, the splitting ([l9|) is exponentially small 
in the system size (at least as small as the estimates (|l4|)- 
( |l6| ) or even smaller). The expression for the splitting 
may be written as 



exp 



(20) 



where R[ n and i? ut are the optimal tunneling trajecto- 
ries from the link (ij) to the inner and outer boundaries, 
and £;„ and £ ut are the corresponding correlation lengths 
(Fig. ^|a). If the system has a "weak spot" — a preferred 
vison tunneling trajectory, the correlations in the imme- 
diate vicinity of that trajectory are most sensitive to the 
choice of the topological sector. Correlations away from 
the optimal trajectory have a weaker splitting between 



the sectors. On the torus, there is no "optimal tunneling 
trajectory" because of the translational symmetry, and 
the correlation- function splitting is given by (p^|)-(^6|), 
independently of the position of the involved operator. 

The above results allow us to estimate the energy split- 
ting in the perturbed Rokhsar-Kivelson model to the first 
order of the perturbation theory. Physically, two types of 
perturbations may be of interest. First, we consider the 
case v/t < 1 in the Hamiltonian (|l|), in which case the 
ground-state energy is no longer zero. The ground-state 
energy may be estimated to the linear order in (1— v/t) as 
(v — t)(H v )nK, where H v is the potential (proportional to 
v) term in the Hamiltonian (Q) , and the average is taken 
in the unperturbed system. H v is a four-point (dimer- 
dimer) correlation function, and from our previous dis- 
cussion it follows that the splitting of (H v )rk in different 
topological sectors is exponentially small in the system 
size with the exponent estimated by (|l4|)-([l6|). Second, 
similar exponential smallness may also be derived for the 
energy splitting under the influence of disorder (again, 
to the linear order in the disorder potential). Note that 
the energy splitting is more sensitive to disorder in the 
vicinity of the optimal vison tunneling trajectory. 

Within our method, the exponentially small response 
of the energy splitting to the perturbation can be de- 
rived only to the lowest order of the perturbation theory. 
Extending those results to higher orders requires infor- 
mation on the excited states at the RK point. This goes 
beyond the scope of this paper, since the excited states 
are not accessible by the Pfaffian technique. 



V. VISON OPERATORS AND CORRELATIONS 

The concept of a vison can be made more transparent 
by introducing a "vison operator" . Given two elemen- 
tary plaquets (triangles) of the lattice and a contour V 
connecting them, we can define the product of two vison 
operators V\ V% as 



V 1 V 2 = (-l) N -=]J(l-2n ij ), 



(21) 



where Nr is the number of dimers intersecting the con- 
tour r, and the last product is taken over all links (ij) 
intersecting T (Fig. ||) (l4| . It can be easily seen that 
deforming the contour T while keeping the end points 
fixed may only change the sign of this operator (depend- 
ing on whether between the old and the new contour lies 
odd or even number of sites). This allows us to define a 
single vison operator (up to a sign; more precisely, the 
vison operator is defined on the frustrated dual lattice, 
with reversing its sign when moved around one site of the 
original lattice). For a single vison operator, the contour 
r is drawn either to infinity (for an infinite system) or 
to the boundary. In a finite system without boundaries 
(e.g. torus), vison operators may appear only in pairs. 

In the Majorana fermion technique (Q), the vison op- 
erator is equivalent to changing the boundary conditions 
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(from periodic to anti-periodic and vice versa, i.e., plac- 
ing a Z2 twist) across the contour T. Thus switching 
topological sectors on a cylinder as described above is 
equivalent to placing a pair of visons in the cylinder open- 
ings. On a torus, switching topological sectors is achieved 
by creating a vison-vison pair, subsequently moving one 
of the two visons around the torus and annihilating it 
back with the other vison. 

Applying a vison operator as defined above to the 
ground state does not produce an eigenstate of the Hamil- 
tonian ([!]). However, preliminary studies show that 
the state thus obtained is close to the actual low-lying 
excitations [^0). A localized excitation (a wavepacket) 
may be classified as either vison-like or non-vison-like de- 
pending on whether it originates a cut T across which the 
boundary conditions are changed. Based on the avail- 
able numerical studies [^0|, we believe that the lowest 
excitations belong to the vison-like sector; this also con- 
forms to the existing field-theoretical description of the 
RVB liquid Q (studies of related models with an exact 
construction of the vison excitations appeared recently 
in Refs. 0, |l0|). We expect that the leading exponential 
asymptotics of long-distance correlators of vison-like op- 
erators is dominated by the topological effect of changing 
the boundary conditions across the cut T, and only finite 
prefactors depend on the details of the involved opera- 
tors. In other words, the correlation length of vison-like 
operators is universal and may be determined, for ex- 
ample, from the pairwise correlation function (V1V2) as 
defined in (^T|). 

We can prove that the vison correlation length equals 
that for Majorana fermions. For the same reason as for a 
hole in the plane (see discussion above) , a vison on an in- 
finite plane produces a zero mode of the matrix A^ . Thus 
the vison correlation length is determined by the decay of 
this zero mode at large distances. The matrix Aij in the 

presence of the vison may be written as Aij — A^' +5Aij , 

where A^ is the matrix of hopping amplitudes without 
the vison and 5 Aij is localized at the cut drawn from the 
vison. Note that the zero mode does not depend on the 
trajectory of the cut (up to a gauge transformation). The 
equation on the zero mode (A 1 - ' 1 + 6A)^o = may be 





FIG. 9: The variational Monte Carlo results of the vison- 
vison correlation function (^l|) in the A direction. The plot is 
in the linear-logarithmic scale, with the distance measured in 
the lattice spacings. The solid dots are the numerical results 
(the circles around dots indicate negative values) . The magni- 
tude of the numerical data should be compared to the analytic 
expression RT 1 ^ exp(— Re R/^a) (empty squares), where £a 
is given by The calculation was performed on the 50x50 
torus with averaging over 10 6 dimer configurations (error bars 
are shown). 



rewritten as ^0 = — [^4^ - ) ] _1 (5A5'o. From this expression, 
it follows that the exponential decay of ^0 away from the 
vison has the asymptotics of LA^] -1 which is exactly the 
Green's functions of Majorana fermions (^|). This fin- 
ishes the proof. Note that this result is consistent with 
our findings for the splitting between topological sectors 
and with their interpretation in terms of vison tunnel- 
ing amplitudes. Furthermore, our result about the decay 
of the zero mode at the point-like vison is a particular 
case of a more general statement about the zero modes 
at even-size holes (see our discussion in Section HI): the 
point-like vison may be considered as a hole of size zero; 
the proof may be extended in a straightforward way to 
cover the general case of an even-size hole. 



FIG. 8: The product of the two vison operators. The visons 
are placed in the shaded triangles. The contour V connecting 
them is shown as the dashed line. The bold links are those 
involved in the operator (|2lj). 



To illustrate our result on the equality between the 
correlation lengths for visons and for Majorana fermions, 
we have calculated numerically the vison-vison correla- 
tion function in the A direction. The calculation was 
performed by the variational Monte Carlo method: the 
vison correlation function was computed as the average 
over a suitable random walk in the space of all dimer 
configurations {jj . Figure || shows the numerical results 
of the calculation (for distances up to 9 lattice spacings, 
on the lattice of 50x50 torus), together with the analytic 
expression R~ 1 / 2 exp(— Re R/^a)- We have no analytic 
expression for the phase of the oscillating part of the cor- 
relation function (arising from the imaginary part of £7*) 
and therefore compare only the overall magnitude of the 
correlation function. The agreement in the magnitude 
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and the change of sign in the vison correlator (a signa- 
ture of the oscillatory behavior) are consistent with our 
analytical prediction of the correlation length (||) . 

To summarize, we have shown that in the dimer model 
(|l|) on the triangular lattice, at the RK point t = v, 
both dimer and vison correlations decay exponentially 
with the distance. Dimer correlations may be expressed 
as finite products of the Green's functions of auxiliary 
Majorana fermions for which the correlation length is 
computed explicitly. We have further shown that the vi- 
son correlation length coincides with that of Majorana 
fermions. The same correlation length also governs the 
splitting between correlation functions in different topo- 
logical sectors (|l4|)-(|l(|). Thus the Rokhsar-Kivelson 
dimer model provides an example of the system where 



the topological ground-state properties of the RVB liq- 
uid may be explicitly verified. 
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P. S. At the final stage of the manuscript prepara- 
tion, we have learned about the recent work |26| , where 
the same dimer model is studied. Many of our results 
in Section [n] overlap with those derived in that work. 
Among other results, Ref. ^6] also contains a more de- 
tailed analysis of the Majorana Green's functions. 
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